#coding:utf-8

import random
import numpy as np
import pylab as p

def teste(a,b,c,x,y,z,r1,r2):
	return pow((x-a),2)+pow((y-b),2)+pow((z-c),2)-pow((r1+r2),2)

def gen_dir(nvar):
	dirs = np.zeros( (nvar+1,nvar) )
	for i in range(0,nvar):
		dirs[i,i] = 1
		
	dirs[nvar,:] = [-1]*nvar
	
	return dirs	
	
def gen_esferas(n,max_r,max_pt):
	xi = [ max_pt*random.random() for j in range(0,n) ]
	yi = [ max_pt*random.random() for j in range(0,n) ]
	zi = [ max_pt*random.random() for j in range(0,n) ]
	
	ri = [ int( max_r*random.random() ) for j in range(0,n) ]
	
	return (xi,yi,zi,ri)

n=5

xi = [2,2,2,2,2]
yi = [2,2,2,2,2]
zi = [2,2,2,2,2]
ri = [6, 7, 2, 3, 4]

'''	
res = gen_esferas(n)
xi = res[0]
yi = res[1]
zi = res[2]
ri = res[3]
'''

FATOR=1.2
MAX_ITER = 6000

def soma_v(v1,v2):
	res = [0]*len(v1)
	for i in range(0,len(v1)):
		res[i] = v1[i]+v2[i]
	return res
	
def interno_v(v1,v2):
	res = 0
	for i in range(0,len(v1)):
		res += v1[i]*v2[i]
		
	return res	
	
def escalar_v(v1,k):
	res = [0]*len(v1)
	for i in range(0,len(v1)):
		res[i] = v1[i]*k
		
	return res	
	
def norma(v):
	res = 0
	for i in range(0,len(v)):
		res+=pow(v[i],2)
	res = np.sqrt(res)
	return res

	
def Ei(i,v):
	ind_i=3*i
	res = 0
	for k in range(0,n):
		ind_k=3*k			
		
		if i==k:
			continue
			
		parc = pow( (v[ind_i]-v[ind_k]),2 ) + pow( (v[ind_i+1]-v[ind_k+1]),2 ) + pow( (v[ind_i+2]-v[ind_k+2]),2 ) - pow( (ri[i]+ri[k]),2 )
		res+=pow(min(parc,0),2)
		
	return res
	
def volume(v):
	max = [-100,-100,-100]
	min = [100,100,100]	
	
	for i in range(0,n):
		ind=3*i
		if v[ind]+ri[i] > max[0]:
			max[0] = v[ind]+ri[i]
		if v[ind+1]+ri[i] > max[1]:
			max[1] = v[ind+1]+ri[i]
		if v[ind+2]+ri[i] > max[2]:
			max[2] = v[ind+2]+ri[i]			
			
		if v[ind]-ri[i] < min[0]:
			min[0] = v[ind]-ri[i]
		if v[ind+1]-ri[i] < min[1]:
			min[1] = v[ind+1]-ri[i]
		if v[ind+2]-ri[i] < min[2]:
			min[2] = v[ind+2]-ri[i]
		
	print abs(max[0]-min[0]),abs(max[1]-min[1]),abs(max[2]-min[2])
	return ( ( ( abs(max[0]-min[0]) )*( abs(max[1]-min[1]) )*( abs(max[2]-min[2]) ) ),min,max )	
	
def Fi(i,v):
	ind =3*i
	res = FATOR*Ei(i,v)
	
	return res
	
def F(v):
	res = 0
	
	max = [-100,-100,-100]
	min = [100,100,100]	
	for i in range(0,n):
		res+=Fi(i,v)
		
		ind = 3*i
		if v[ind]+ri[i] > max[0]:
			max[0] = v[ind]+ri[i]
		if v[ind+1]+ri[i] > max[1]:
			max[1] = v[ind+1]+ri[i]
		if v[ind+2]+ri[i] > max[2]:
			max[2] = v[ind+2]+ri[i]			
			
		if v[ind]-ri[i] < min[0]:
			min[0] = v[ind]-ri[i]
		if v[ind+1]-ri[i] < min[1]:
			min[1] = v[ind+1]-ri[i]
		if v[ind+2]-ri[i] < min[2]:
			min[2] = v[ind+2]-ri[i]
		
	res += ( ( abs(max[0]-min[0]) )*( abs(max[1]-min[1]) )*( abs(max[2]-min[2]) ) )
		
	return res
		
def main():
	delta = 5.0
	epson = 1e-6
	it = 0
	
	nvar = 3*n
	dirs = gen_dir(nvar)
	x0 = [0]*nvar
	
	f = F(x0)
	
	while (delta>epson and it<=MAX_ITER):
		print it,f
		found_dir = False
		for d in dirs:
			xk = soma_v(x0,escalar_v(d,delta))
			f_ = F( xk )		
			
			if f_ < f:
				x0 = xk
				f=f_
				found_dir=True				
				break		
				
		if not(found_dir):
			delta/=1.2
			
		it+=1
	
	print it,delta
	return x0
	
	
x = main()


#Print Scilab plot

print "v=[",
for i in range(0,len(x)):
	if i%3==0 and i!=0:
		print ";",
	print x[i],
print "]"

print "ri=",ri

resp = volume(x)
print "caixa=[ ",resp[1],"],[",resp[2],"]"
print "Volume=",resp[0]
